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' Abstract 

1^ Gamma-ray burst analyses at neutrino telescopes are typically based on diffuse or stacked 

'""' (z. e., aggregated) neutrino ffuxes, because the number of events expected from a single burst 

cn is small. The interpretation of aggregated flux limits implies new systematics not present 
for a single burst, such as by the integration over parameter distributions (diffuse fluxes), or 



00 by the low statistics in small burst samples (stacked fluxes). We simulate parameter distri- 

butions with a Monte Carlo method computing the spectra burst by burst, as compared to 
a conventional Monte Carlo integration. With this approach, we can predict the behavior 
(^ of the flux in the diffuse limit as well as in low statistics stacking samples, such as used in 

^~~^ recent IceCube data analyses. We also include the flavor composition at the detector (ra- 

tio between muon tracks and cascades) into our considerations. We demonstrate that the 
. !^ spectral features, such as a characteristic multi-peak structure coming from photohadronic 

/\ interactions, flavor mixing, and magnetic field effects, are typically present even in diffuse 

C^ neutrino fluxes if only the redshift distribution of the sources is considered, with ^ ~ 1 dom- 

inating the neutrino flux. On the other hand, we show that variations of the Lorentz boost 
can only be interpreted in a mo del- dependent way, and can be used as a model discrim- 
inator. For example, we illustrate that the observation of spectral features in aggregated 
fluxes will disfavor the commonly used assumption that bursts with small Lorentz factors 
dominate the neutrino flux, whereas it will be consistent with the hypothesis that the bursts 
have similar properties in the comoving frame. 
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1 Introduction 

Neutrino telescopes, such as IceCube |1| or ANTARES [2], are expected to reveal the nature 
of the sources of the highest-energetic cosmic rays. If a substantial amount of protons (or 
heavier nuclei) are accelerated in the sources, charged pion production is expected via pp and 
P7 interactions, and therefore a neutrino flux. There are numerous candidate sources, see 
Ref. [3] for a review and Ref. El for the general theory. We focus on the prompt emission of 
gamma-ray bursts (GRBs) in this study, where photohadronic (^7) interactions are expected 
to lead to a significant flux of neutrinos |5|. If GRBs are the sources of the highest-energetic 
cosmic rays, the expected neutrino flux is just around the corner. In fact, recent IceCube-40 
data (referring to the 40 string configuration of IceCube) start to test this paradigm |6l[7] , 



based on analytical estimates of the neutrino fluxes, see Refs. pip 10 and Ref. 11 for 
the application to IceCube data. On the other hand, numerical calculations indicate that 
additional processes in the photohadronic interactions, flavor mixing, and magnetic field 



effects affect normalization and shape of the expected neutrino fluxes, see Refs. 12 14 . In 



Ref. 14 , we have made the impact of these effects on the original Waxman-Bahcall (WB) 
flux shape [5] very explicit. For example, a characteristic multi-peak structure from high- 
energy processes in the photohadronic interactions, magnetic field effects, and flavor mixing 
has been predicted for the same astrophysical assumptions as for the WB flux. 



While Ref. 14 is essentially based on one set of parameters only, we focus in this study on 
aggregated neutrino fluxes. Because of the low statistics expected from a single burst, any 
state-of-the-art analysis is based on such an aggregated flux, which can either be a diffuse 
flux, for which individual sources are not resolved, or a stacked flux, for which the neutrino 
fluxes from individual sources are expected to be correlated with their gamma-ray counter- 
parts. The interpretation of aggregated flux limits, of course, implies new systematics not 
present for a single burst, such as by the integration over parameter distributions (for diffuse 
fluxes), or the assumptions for unknown burst parameters (stacked fluxes). In the stacking 
case, the time- and directional-wise correlation can be used to suppress backgrounds, and 
the fluence of the gamma-rays can even be used to compute the expected shape of the 
neutrino spectrum [61. However, since for many bursts the parameters needed for such an 
estimate are unknown, such as the redshift or Lorentz factor, "standard values" are often 
used. It is one of the motivations for this study to critically review these standard values 
from a theoretical perspective, in order to reduce the systematical errors in the interpreta- 
tion of future analyses. For instance, depending on the method, a too large standard value 
for the redshift may easily lead to an overestimated neutrino flux. In addition, we discuss if 
the spectral multi-peak features are present in aggregated fluxes. Note that experiments as 
IceCube have limits with optimal sensitivities in specific energy ranges, see e.g. discussion 



in Refs. 15, 16 . Therefore, it is relevant if these spectral features are present in aggregated 



fluxes, especially in a diffuse flux. 



Predictions for diffuse neutrino fluxes have, for instance, been made in Refs. 12,17,18 . For 
a diffuse flux, typically at least the redshift distribution is integrated over, which is assumed 
to follow the star formation rate in some form. Apart from the redshift, interesting parame- 
ters which vary from burst to burst are Lorentz factor, luminosity, gamma-ray fluence, and 



variability timescale. For example, in Ref. 18 , relying on conventional fireball phenomenol- 



ogy, a log normal distribution for the Lorentz factor is assumed. In the conventional fireball 
model, the main contribution to the neutrino flux is assumed to come from small Lorentz 



factors, see Refs. 19,20 for discussions on the parameter space constraints. This conclusion 
is based on the estimate of the size of the interaction region from the variability timescale, 
which strongly depends on the Lorentz factor. However, in a modified version of the fireball 
model including with the hypothesis that the bursts have similar properties in the comoving 
frame (shock rest frame) the dominance of the small Lorentz factors is no longer given, as we 



will demonstrate. In fact, a recent analysis by Ghirlanda et al. 21 , which appeared during 



completion of this work, favors this hypothesis from the simple argument that the variations 



of the bulk Lorentz factor imply the Amati 22 and Yonetoku 23 relations, while the bursts 



are expected to be similar in the comoving frame. In order to clarify the model-dependence 
(c/.. Sec. [2]), we start with the minimal ingredients (and input parameters) to neutrino pro- 
duction needed to numerically reproduce the analytical estimates in Refs. [SllGp] . We clearly 



separate these from model-dependent assumptions (Sec. 2.2), which has the advantage that 
one can clearly see where the model-dependent hj^otheses enter. The main questions which 
we pursue in this context are: Assuming theoretical expectations for parameter distribu- 
tions, from which parameters do the main contributions to the neutrino flux come from? As 
a consequence, what value should one assign to a burst with unknown redshift or Lorentz 
factor, such as in a stacking analysis? What are the assumptions going into that conclusion, 
in particular, which assumptions hold in a model-independent way? In order to discuss 
these questions in a coherent and transparent way, we discuss each of the relevant input 
parameters separately in Sec. [3j 

In principle, a neutrino telescope cannot only measure the neutrino flux, but also distinguish 
different flavors by different event topologies |24]. In the simplest case, one may use two 
topologies: tracks (mainly from z/^) and cascades (mainly induced by v^ and z/^) to construct 
an observable muon track to cascade flavor ratio as a flavor dependent quantity. Note that 
this flavor ratio is defined in the spirit of cascade searches from extragalactic neutrino 



sources, which has been very recently conducted in Ref. 25 for IceCube-22. From the 
charged pion vr — )■ /i — )• e decay chain, the flavor composition at the source is given by 
z/e : t'u : z/t- = 1 : 2 : ("pion beam source"), and, for GRBs, changes to : 1 : at higher 



energies 13,26] ("muon damped source") because of the synchrotron cooling and decay of 
the intermediate muons, pions, and kaons. Although one may not expect to observe this 
effect with high statistics on the timescale of ten years of full IceCube running due to the 
lower effective area for cascades, some conclusions on the astrophysical source (such as the 



magnetic field |27]) or the particle physics during neutrino propagation (see, e.g., Ref. 28 
for a short review) may be obtained at IceCube or future potential upgrades. Since the flux 
normalization drops out of this flavor ratio, it ought to be relatively robust with respect to 



astrophysical uncertainties 29 . In this study, we do not only discuss the muon neutrino 
flux at the detector, but also the flavor ratio. In particular, we will demonstrate that the 
aggregated neutrino flavor ratio, though perhaps measurable with much lower statistics, is 
a more reliable prediction for GRBs than the flux itself; see Sec. [3] 

Finally, we discuss the systematical error on the quasi-diffuse flux introduced by the stack- 
ing statistics, i.e.., the number of stacked bursts, in Sec. |4} Note that this systematical 
error is important for the interpretation with respect to the paradigm that GRBs are the 



sources of the highest energetic cosmic rays. It is also illuminating to show what the main 
qualitative differences between stacked and diffuse fluxes are, and where we expect the bor- 
derline. Throughout this paper, we use a Monte Carlo method to compute the neutrino 
fluxes. Instead of Monte Carlo integration over the parameter distributions, we compute the 
neutrino spectrum burst-by-burst, sampling over the input parameter distributions. This 
method has the advantage that it converges to conventional (Monte Carlo) integration in 
the limit of a large number of bursts, such as (9(10 000) bursts expected in the observable 
universe over a timescale of ten years. As we shall see, this limit then also represents a 
diffuse flux. On the other hand, small sample sizes can be chosen to simulate the impact of 
stacking analyses of e.g. 0(100) bursts, as used in Ref. [6]. 

2 The source model 



Here, we describe the model we use for neutrino production, see also Ref. 14 . First, in 



Sec. |2.1[ we give a short summary of the model and discuss the translation into observables. 
Details can be found in App. |X| where also a reference parameter set reproducing the 
original WB shape is derived, which we use in the following sections, and the effects of the 
photohadronic high-energy processes and secondary cooling are shown quantitatively for this 
example. We also discuss the limitations of the model in greater detail there. This model 
description is kept as independent as possible from any underlying astrophysical gamma-ray 
burst model. In the next step, in Sec. |2.2[ we discuss the connection with observables and 



we demonstrate where model-dependent relationships enter. 
2.1 Short model summary 



The idea of this model is to describe the analytical approaches in Refs. (5lp, 11 as closely 
as possible, which are used for state-of-the-art IceCube stacking analyses, while introducing 
more realistic descriptions for energy losses and escape required to describe magnetic field 
and flavor effects. As indicated above, we clearly separate the (minimal) set of ingredients 



needed for the description of the neutrino spectra (Sec. 2.1) and the methods to estimate 



these quantities in a model-dependent way (Sec. 2.2). The ingredients of our model closely 



follow Ref. 27 , where the main differences is the origin of the target photons (which come 
from the synchrotron radiation of co-accelerated electrons/positrons in Ref. [27]). Note 
that in this work, we do not revise the absolute normalization of the expected neutrino flux, 
as it may be expected from the gamma-ray or cosmic ray counterpart, but focus on the 
theoretically anticipated relative contributions from different bursts. Instead we use for the 
expected normalization of the (diffuse) flux [30| (updated in Ref. fsTl) 

El(j)^ = 0.45 ■ 10"^ -^ GeV cm'^ s"^ sr-^ (1) 

per neutrino species [ue, i^fi, or z/^ from the pion decay chain) before flavor mixing, or, to a 
first approximation, for any flavor after flavor mixing. This normalization depends on the 
fraction of the proton energy going into pion production /^, for which we assume /,r — 0.2. 



We will examine the assumptions going into this frequently used estimate elsewhere 32 . 



The photon spectrum is assumed to fit the observed spectral energy distributions of GRBs 
described as a broken power law with a spectral break, parameterized in the shock rest 
frame (SRF)/comoving frame by 



7, break / 



pi < f' < f' 

^7,min — ^ ^ ^7, break 



^7(^0 = C'7 ■ S f —^ V ' f' <f' < f' 

' ' I e' V , / ^7,break -i ^ ^ fc7,max 

\ 7, break / 

else 

where C is a normalization factor in GeV~^ cm~'^. Typical values are a^ ~ —1 and /3^ ~ —2. 
The proton spectrum is assumed to be a cut-off power law with a spectral index expected 
from Fermi shock acceleration and a cut-off energy constrained by the synchrotron losses of 
the protons: 




N'^{E'^) = C'^- I V^eV; ^^y\ S;,%.J -p - -p,min ^ ^3) 



where we typically choose Up = —2 and C is the normalization of the spectrum (in units 
of GeV-^cm-3). 

Photohadronic interactions lead to charged pion, kaon, and neutron production. These 
particles decay further through weak decays into muons and neutrinos. Our photohadronic 
interaction model is based on Ref. [33] (model Sim-B), based on the physics of SOPHIA |34] , 
including direct production, different resonances, and high energy processes. The helicity 
dependent muon decays are taken into account as in Ref. p!3|. We assume that the source 
is optically thin to neutrons, which means that neutrons produced by the photohadronic 
interactions will decay and lead to an additional electron antineutrino flux. For all produced 
charged secondaries, we consider synchrotron energy losses and escape via decay, which leads 
to a loss steepening of the spectrum at 



depending on the properties of the secondaries (mass m and rest frame lifetime tq) and the 
value of B'. These spectral breaks, together with high energy processes in photohadronic 
interactions, translate into a characteristic multi-peak structure in the neutrino spectra, see 



App. A. 2 Although not all features of the sources may be described by our method (see. 



e.g., Ref. |T^), it has the advantage that the results can be directly related to Refs. [5||8||TI]. 

To summarize, we show in Fig. [T] a flowchart describing the computation of the neutrino 
spectra in the SRF. Here, Q'{E') (rectangular boxes) stand for (injection) spectra per time 
frame and N'{E') (boxes with rounded corners) for steady spectra derived from the balance 
between injection and losses or escape. Dashed arrows stand for solving the steady state 
differential equation balancing energy losses and escape with injection. Weak decays and 
photohadronic interactions are denoted by the horizontal lines. The resulting neutrino 
spectra are marked by a double border and a colored background. Below the boxes we 
indicate by which parent the neutrinos are produced. 
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Figure 1: Flowchart describing the model (in SRF). The functions Q'{E) denote (injection) spectra 
per time frame [(GeV cm^ s) ] and N'{E) steady spectra [(GeV cm^) ] derived from the balance between 



injection and losses or escape. Dashed arrows stand for solving the steady state differential equation Eq. (32 1 



The transformation of tlie injection spectrum of the neutrinos Q'^,^ (in units of GeV~^ cm~'^ s~^) 
from a single source into a point source or quasi-diffuse flux (in units of GeV~^ cm~^ s~^ [sr~^]) 
at the detector 0^ before flavor mixing is given by 



A^ 



1 



Audi 



Q'u, 



E, 



— El. 



(5) 



Here, iV is a (model-dependent) normahzation factor depending on the volume of the in- 
teraction region and its possible Lorentz boost with respect to the observer, which will be 

[z) is the luminosity distance of the 



discussed later in Sec. 2.2, Here, ^^(z 



:i + z) 4 



source computed in standard (fiat) FLRW cosmology with the values from Ref. 35 and 
'^com(^) is the comoving distance jj From Eq. ([s]), one can read off that the neutrino flux 
E'^(f) oc l/c?| independent of the model if only the redshift dependence is considered. In 
order to compute the 0/3 (/3 = e, /i, r) neutrino flux at the detector, we sum over all these 
initial neutrino fluxes of flavor a weighted by the usual flavor mixing 



Ei^" 



\u, 



(Si\ 



(6) 



We use sin^ 9i2 = 0.318, sin^ 6*23 = 0.5, and sin^ ^13 = for the sake of simplicity (see. 



e.g., Ref. 36 ; 6*13 = is compatible with recent T2K hint at 2.5a |37)), leading to flavor 



^Note that the flux, as we define it per energy, time, and area, scales ex l/dcomj where as the energy flux 
density, defined as energy per time frame and area, scales oc l/d|. 



Parameter Units 



Description 



Typical 



Ref. 



B' 

r 

z 



G (Gauss) Magnetic field strength 
1 Lorentz boost factor 

1 Redshift 



lO^G.-.lO^G 
100 .. . 1400 
0.02... 6 



300 kG 

102.5 

2 



'7, break 



a. 

13, 



a. 



7 



keV 
1 
1 
1 



Break in photon spectrum (SRF) 0.1 keV, 

Photon index below break —1 

Photon index above break —2 

Proton spectral index —2 



..30keV 14.8 keV 
-1 
-2 
-2 



7] 

N 



1 
a.u. 



Acceleration efficiency 
Normalization/luminosity 



0.1 

arbitrary units 



0.1 
Eg. ^ 



Table 1: Parameters of our model, units, description, typical values used in this study, and a reference 



parameter set (chosen in Sec. A. 2) 



equipartition between z/^ and Vr at the detector. We sum over neutrinos and antineutrinos, 
i.e., if we refer to "z/, 



^l 1 



we mean the sum of the t'^ and v^ fluxes. 



Apart from the flux, the simplest observable may be a flavor ratio in which the detector 
properties are taken into account. Except from muon tracks, the event topologies currently 
discussed in IceCube are cascades 25 . If we assume that electromagnetic (from z/g) and 
hadronic (from Vt) cascades do not need to be distinguished, a useful observable is the ratio 
of muon tracks to cascades 



38 



R 



+ 



(7) 



Note that neutral current events will also produce cascades and u^. will also produce muon 
tracks in 17% of all cases, which, in practice, have to be included as backgrounds. In 
Ref. 25 , the contribution of the different flavors to the cascade rate for a i?~^ extragalactic 



test flux with equal contributions of all flavors at the Earth was given as: electron neutrinos 
40%, tau neutrinos 45%, and muon neutrinos 15% (after all cuts). This implies that charged 
current showers dominate and that electron and tau neutrinos are detected with comparable 
efficiencies, i.e., that Eq. ([7]) is a good first approximation to discuss flavor at a neutrino 
telescope. The benefit of this fiavor ratio is that the normalization of the source drops out. 
In addition, it represents the experimental flavor measurement with the simplest possible 
assumptions. For a pion beam source, where at the source z/g : t'/x : i^r = 1 : 2 : 0, at the 
detector R is 0.5, for a muon damped source, where at the source z/g : z/^ : z^,- = : 1 : 0, 
at the detector R is 0.64, and for a neutron beam source (neutrinos from neutron decays), 
where at the source z/g : z/^ : z/^ = 1 : : 0, at the detector R is 0.28. 

Our main model parameters are summarized in Table [ij In this table, the set of parameters 
in the last column "Ref." has been chosen as one possibility to reproduce the original WB 



shape in the analytical version, see App. A. 2 In the following, we use Ercj)^ from Eq. m 



representative for the muon tracks in the neutrino telescope, and R in Eq. (M), representative 
for the muon track to shower ration in the neutrino telescope, as observables, i.e., we will 
show the results for these observables where apphcable. 



2.2 Model dependence and connection to gamma-ray observations 

In order to demonstrate how burst to burst variations are to be interpreted in terms of dif- 
ferent models, we discuss three different approaches: cannonball model-like phenomenology, 
a general A^ zone fireball model with properties similar in the SRF, and conventional fireball 
phenomenology. In the latter two cases, we illustrate how the input parameters needed (c/.. 
Table [I]) can be related to gamma-ray observations, in particular, how the photon spec- 
trum normalization C in Eq. (|2| and the proton spectrum normalization C in Eq. (JSJ) are 
obtained explicitely. 

Cannonball-like phenomenology (CB) 

Suppose that the acceleration region emits isotropically in the bulk rest frame, and this 
whole region is ejected by the engine with a Lorentz boost F, such as a cannonball (CB), 



see Ref. 39 for a specific modeln Then the isotropic emission region will be boosted into a 
solid angle oc l/F^ by simple relativistic dynamics, and the normalization factor in Eq. dsl) 
can be written as 

N = V'V^ = —B!^V\ (8) 

As a consequence, Elcf) oc F^ if the cannonballs are otherwise alike in terms of the photon 
and proton spectra in the bulk rest frame, i.e., the Lorentz boost F is entirely a feature of 
the engine. This is, of course, a very specific assumption, which does not take into account 
the origin of the target photon field, as in Ref. p9|jj Finally, note that the probability 
to observe such a burst will be oc 1/F^, which is determined by the probability that the 
observer is within the solid angle of the jet. 

Fireball model with bursts alike in shock rest frame (FB-S) 

If the gamma-ray burst comes from a relativistically expanding fireball (FB), we can assume 
that the acceleration region expands isotropically in the source (engine) frame. This is even 
a good approximation if there is a jet with the opening angle 9 > 1/F, since the Lorentz 
boost will "overwrite" the beaming by the Lorentz transformation into the observer's frame 
as long as the observer is close enough to the jet axis. Therefore, the interaction region 
is typically assumed to be isotropic in the source frame, and the beaming does not have 
an effect on the quantities we are interested in. On the other hand, it is implied that the 
probability to observe such a GRB is proportional to the solid angle of the jet ^^/2, which 
may be correlated with F. 

The normalization factor in Eq. ([s]) cannot simply be obtained as in Eq. (Is]), because the 
emission region is not spherically symmetric in the SRF, but in the source frame. Since it is 

^Note that, in general, a Doppler factor depending on the viewing angle has to be used, see e.g. Ref. |4]. 



'In fact, in Ref. 39 a thermal photon field which is isotropic in the source frame is assumed. Therefore, 
the photon energies in the SRF are proportional to T (and, in addition, the Lorentz boost and Doppler 
factor have to be distinguished). Here, for the sake of simplicity, we assume that the target photon field is 
isotropic in the bulk, as for synchrotron photons, but this does not affect our argument. 



a frequently used approach to estimate the neutrino spectrum from the gamma-ray signal, 
we demonstrate how this normalization is typically performed. Because our model does 
not describe the neutrino production in a time-resolved way, it makes sense to relate the 
neutrino production to the gamma-ray fluence during the burst. The quantity of interest is 
the bolometric fluence St,oi (in units of ergcm"^). In practice, of course, only information 
from specific wavelength bands may be available. It can be used to calculate the isotropic 
bolometric equivalent energy -Eiso,boi (in erg) in the source (engine) frame as 

-E'iso,boi = jz — ; — r 5'boi • (9) 

(1 + ^) 

This energy corresponds to the total emitted photon energy at the source if the emission is 
isotropicn It can easily be boosted into the SRF by -Ej'so.boi — -^iso,boi/r- Assuming that the 
emitted photons are coming from synchrotron emission of electrons (or mainly interact with 
electrons), the amount of energy in electrons and photons should be roughly equivalent. 
Therefore, if the photons carry a fraction e^ of the total energy -Eiso,tot5 we have 

-E'iso.tot = Cg ■ -^iso,bol • (10) 

In order to compute the photon and proton densities in the SRF, it turns out to be useful to 
define an "isotropic volume" Viso, which is the volume of the interaction region in the source 
frame assuming isotropic emission of the engine. Similarly, V^^^ is an equivalent volume in 
the SRF where only the radial direction is boosted. Because of the intermittent nature of 
GRBs, the total fluence is assumed to be coming from N such interaction regions. Now one 
can determine the normalization of the photon spectrum in Eq. ^ from 

E' 

I i\tI f '\ J ' iso.bol /I 1 \ 

eN^{e)de = ^^tttt ■ (H) 

iso 

if one assumes that the target photons actually escape from the source on the timescale 
considered. Similarly, one can compute the normalization of the proton spectrum in Eq. (Is]) 
by 

where /e is the ratio between energy in electrons and protons (/7^: baryonic loading). Note 
that in the end, we will obtain the neutrino flux per time frame per interaction region 
from Eq. ([5| with N = V(^^. Given that the magnetic field carries a fraction eb of -Eiso,tot, 
one has in addition 



ttI ^B iso,bol 7-)/ /o ^B isOjbol /t o\ 

Typical values used in the literature are /e ~ eg ~ es — 0.1 (see, e.g., Ref. [TT]). 



^In fact, for this computation, we do not need the actual (coUimation corrected) emitted energy, which 
depends on the beaming of the jet, because any beaming effect will cancel out in our computation of the 
proton and photon densities. 



From Eqs. (11) and (12), together with Eq. ([5]), one can easily read off that if different 
bursts have the same characteristics in the SRF, such as A^, V^^^, and E[^^^^i, there is no 
dependence of the neutrino flux on F, and El(j) oc F^. As for the cannonball above, 
this case may be plausible if the properties in the bulk are alike, but the Lorentz boost is 
determined by the engine. This hypothesis is consistent with Ref. (2l], assuming that the 
bursts are similar in the comoving frame. We will refer to it as "FB-S" (fireball-shock) in 
the following. 



Conventional fireball phenomenology (FB-D) 

Here, we describe how the unknown parameters in the previous A^ zone model are related 
to conventional fireball phenomenology. The gamma-ray emission is assumed to originate 
from the collision of internal shells of the thickness Ad' at the collision radius Re, see 
Refs. 40 , 41 for reviews. These parameters determine the size of the interaction region, 
i.e.. 



in our language 



v 

' ISO 



AirRl-Ad' 



47rR'c-T-Ad 



(14) 



The time variability of the gamma-ray signal can be described by many such collisions. In 
turn, the variability timescale t^ of the signal at the observer can be used to draw conclusions 
about the geometry, leading to the well known estimates 



Rc^2r^c 



U 



[l + z] 



Arf' ~ F ■ c 



t: 



'l + z) 



Rc 

2F' 



(15) 



where the effects of the cosmological redshift are taken into account. From Eq. (14) we can 
then estimate the size of the interaction region as 



V' ~ 

ISO 



Atc 2 F^ c 



Zy 



[l + Z] 



F-c 



6-); 



'l + z) 



(xT 



(16) 



In addition, if there are N collisions, we can estimate that A^ ~ Tgo/t^, where Tgo is the time 
during which 90% of the total energy is observed. Consequently, the total neutrino fluence 
is obtained by multiplication with t^ N = Tqq. Eq. (16) implies that for fixed ty, the larger 



F, the larger the interaction region, and the smaller the photon density in Eq. (11), which 
directly enters the fraction of fireball proton energy lost into pion production /^ oc F~^ M , 
or, consequently, El(f) oc F~^. Therefore, the main contribution to the neutrino flux is often 



believed to come from bursts with small Lorentz factors, see discussions in Refs. (8,19,20 
As a consequence, recently observed bursts with high F's by Fermi-LAT would, despite 
their high luminosities, not lead to significant neutrino production. We will refer to this 
hypothesis as "FB-D" (fireball-detector) in the following, since the conclusions are based on 
the uncorrelated variation of the detected quantities. 



We can now compute the magnetic field in Eq. (13) using Eq. (15) as 



B' ~ 220 
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in consistency with Refs. [8,10 ^ This means that the reference values chosen in Table ^ in 
particular, the magnetic field, are consistent with the values of Tgg, ty and -Eiso.boi suggested 
by this formula. 

Comparison among approaches 

Comparing the different approaches, CB and FB-S are attractive because variations of the 
bulk Lorentz factor F of otherwise similar bursts in the comoving frame lead to natural 
descriptions of observed empirical relationships, such as the Amati and Yonetoku relation- 



ships, which are obtained from simple relativistic transformations; see, e.g., Refs. 21,42 
On the other hand, FB-D is the most frequently adapted alternative in neutrino physics. 
From the discussion above, it is clear that the conclusion that small Lorentz factors domi- 
nate the aggregated neutrino fiux only holds in FB-D, based on the geometry estimate in 



Eq. (16). However, the exact shell widths are unknown. More importantly, this conclusion 



implies that ty is held fixed and therefore considered as input. However, t^ is a not even 



well-defined observable. For example, from the discussion in Ref. 43 , it is clear that there 
can be several superimposed variability components. In addition, the variability depends 
on the wave length band and the fastest variation may not even be observed because of 
the instrumental cutoff|_] On the other hand, the alternative FB-S corresponds to the very 
recently found correlations in Ref. 21 , which translates into -E^^i/ oc F^. Of course, there 



may be additional interesting correlations, such as F oc tjy for fixed V^^^; cf., Eq. (16). 



3/5 



Since it is not clear which parameters are to be considered as input (the properties at the 
source or at the observer, such as ty), we will test three hypotheses: Elcj)^, oc F^ (CB), 
£"^01, oc F^ (FB-S; parameters in SRF fixed), and E^cp^, oc F~'^ (FB-D; parameters at de- 
tector fixed). The truth may lie somewhere in between, but will require the further test of 
empirical correlations. 

3 Diffuse neutrino fluxes 

Diffuse neutrino fluxes are generated by many bursts, which are assumed to be not resolvable 
individually. On the other hand, stacking analyses are based on a "flnite" number of bursts, 
typically a few hundred bursts observed in gamma-rays. From a stacking limit, a quasi- 
diffuse limit can be computed if the total number of bursts per year is roughly known. This 
limit will obtain a systematical error coming from the statistics of the stacking sample, see 
Sec. |4J Note that there may be also selection effectqjor uncertainties in the total number of 
bursts, which we do not discuss in this work. In practice, it is assumed that (9(1000) GRBs 
per year happen in the observable universe, see, e.g. , Ref. (ol. Therefore, even a diffuse flux 



^In fact, different versions of this formula uses L^aso- However, as long as A^ = Tgo/ty, these results are 
identical apart from the redshift effect on ty, which we take into account. 



^See, e.g., Figs. 2 and 3 in Ref. 43 , where in some cases the a dip seems to be present below the cut-off 
frequency. 

^One example for such a selection effect would be the threshold effect for the minimal detectable flux, 
which is important for the detection in photons. In App. IB] we briefly discuss the consequence of this 
selection effect. 
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is based on a finite number of GRBs, wliicli are, for ten years of operation of the neutrino 
telescope, C(IOOOO) GRBs contributing to the diffuse flux. Therefore, we use the same 
method to compute diffuse and quasi-diffuse fluxes: we generate a probability distribution 
for one parameter, say redshift z, and compute n neutrino fluxes with parameters drawn 
randomly from this distribution. The aggregated neutrino flux is then summed over these 
n individual burst fluxes, and, after that, normalized to our reference flux described by 
Eqs. ([I]) and (36) in terms of the total energy flux of the neutrinos. Note that this method 



takes into account the relative contribution from the individual fluxes, but the absolute 
normalization of the aggregated flux is not an outcome of our approach; it will be discussed 



elsewhere 32 . For the sake of simplicity, we do not distinguish between short and long 
bursts in this analysis, while we choose typical parameters of long bursts. 

The method has two advantages: First of all, diffuse and quasi-diffuse fluxes are computed 
with the same method, which means that, by varying ra, the borderline between diffuse and 
quasi-diffuse limits can be illustrated, and low statistics effects can be quantified, see Sec. |4j 
Second, it may be more realistic than a pure integration of the flux, especially if there are 
observable effects between (9(10 000) and infinitely many bursts contributing to the diffuse 
flux. As we will demonstrate later, n = 10 000, which we use throughout this section, is 
already a very good approximation for the diffuse limit. But this is an outcome of our 
analysis, not an input assumption. The results in this section are also representative for the 
integrated flux from a stacking analysis with very high statistics. Note, however, that the 
neutrino effective area of the neutrino telescope may lead to selection effects of individual 
bursts, which we do not discuss in this work. 

Since each of the parameters, as we will discuss in the following, has different theoretical 
issues to be addressed, we vary only one parameter at the time for clarity, and keep the 
others fixed. We consider the main input parameters for our model, see Table [l} z, B', T, 
and the photon spectral shape. For each parameter, there are mainly two challenges: 

1. The description of the (theoretical) distribution of the parameter. 

2. The interpretation of a different value of the parameter in terms of the relative con- 
tribution to the total flux. 

For example, the redshift z is a relatively simple case: The distribution is typically assumed 
to follow the star formation rate in some form, and the consequences of a different parameter 
value are independent of the source model if the GRBs are otherwise equivalent at the source. 



A very different example is F (see discussion in Sec. 2.2): The theoretical distribution 
of the parameter is not directly known, and the consequences of a different value of F 
depend on the assumptions for the model. Therefore, it is clear that variations of z and F 
cannot be compared at the qualitative level. Before we vary a parameter, we also show the 
consequences of different values explicitely. 

3.1 Redshift 

Here, we assume that all GRBs are alike at the source and the different fluxes and spectra 
are generated solely by the different redshifts of the bursts. We show in Fig. [2] the impact of 
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Figure 2: Neutrino flux (left) and flavor ratio (right) for the reference values from Table fl] where the 
redshift z is varied. The burst with z = 2 is normalized to the WB spectrum, the other bursts scale with 
redshift as given by Eq. ([s]). 

the redshift on the neutrino flux (left) and flavor ratio (right) for the GRB with the reference 
values from Table [l| with z varied. From the discussion after Eq. (|5|, we have -E^^i/ oc d^^, 
i.e., bursts with small redshifts contribute most. Since the energy i^^, oc (1 + z)~^, there is 
hardly any shift of the spectrum, as it can best be seen in the right panel. Nevertheless, it 
is not clear if the contributions from different bursts will lead to some averaging as long as 
the distribution of bursts is not considered. 

Since long GRBs are assumed to originate from the death of massive stars, it is natural 
to assume that the GRB rate traces the star formation rate. For example, Hopkins and 
Beacom use a piecewise steady approach to parameterize the star formation rate density 
(per comoving volume) as [44] : 

r (1 + z)3-44 0<z< 0.97 

p,{z)(x< 10^-09 ■ (1 + z)-o-26 0.97 <z< 4.48 . (18) 

[ lQ^-^^-{l + z)-'^-^ 4.48 < z<6 

As a maximum value for the redshift, Hopkins and Beacom use a cut-off at ^max = 6 because 
of the low statistics at higher z. Furthermore, we use Zmm — 0.02 from the scale the universe 
becomes homogeneous (about 100 Mpc). We assume that the redshift distribution of the 
bursts follows (see, e.g., Ref. |45j) 



dN ^, , , , dV/dz 
— oc£(.)p.(^)y^ 



(19) 



Here, we have used the comoving volume per unit redshift dy/d2(l + z) ^, where the 
factor 1 + z accounts for the cosmic time dilation. In addition, we use the correction 
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Figure 3: Distribution of 10 000 bursts AN jAz as a function of redshift (histograms) and relative con- 
tribution of the individual GRBs d^ dTV/dz (solid curves). The dashed curves show the exact distribution 
functions. The left panel uses the star formation rate from Hopkins and Beacom [44 with the correction 



E[z) from Kistler et al. 45 , the right panel shows the "star formation rate 1" from Porciani & Madau 46 
without correction. 



Ei^z) = Sq{1 + z)^'^ introduced in Ref. 1451 to account for tlie fraction of stars resulting 
in GRBs (strong evolution case). We show the distribution of 10000 bursts in redshift in 
Fig. [sl left panel, for the distribution in Eq. (19) as histogram. For comparison to a very 
different model, the distribution for the model "SFl" by Porciani and Madau [46| (without 
correction factor, i.e., star formation rate evolution) is shown in the right panel of Fig. Isl 

The "typical value" for the redshift of a GRB according to the star formation rate in 
Eq. (18) is roughly around z = 1 — 2. With the correction factor S{z) these values are 
shifted to higher redshifts with a broad peak around z = 1 — 5 (left histogram). In the right 
histogram (without correction factor), the peak is at about z = 1 — 3. However, do the 
GRBs around these values also contribute most to the neutrino flux? For this question, we 
show in both panels of Fig. [sjthe function l/(i| dN/dz measuring the relative contribution 
of individual bursts, which includes both the redshift distribution function dN/dz and the 
weight function l/diizy scaling E^cf) {cf., discussion after Eq. ([s])). Clearly, the weight 
function prefers smaller values of z, whereas the redshift distribution pulls to larger values, 
leading to a local narrow maximum at z = 1. From this generic argument, we conclude 
that the "typical GRB" dominating the neutrino flux has 2; ~ 1, rather than larger values 
of z often used in the literature. If the redshift of a GRB is not measured, it may therefore 
be the best estimate for the neutrino flux computation. For instance, choosing 2; = 2 or 
z = 3 for the unknown redshift may, depending on the method, lead to an overestimate of 
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Figure 4: Neutrino flux (left) and flavor ratio (right) for the redshift distribution according to Hopkins 
and Beacom (HB), as depicted in Fig. l3J left panel. The thick solid curves show the integrated flux, the thin 
solid curves the contributions from individual redshift ranges (left panel only) . The dashed curves show the 
reference values from Tabic [l] The thin dashed curve in the left panel is the WB reference, the total fluxes 
are normalized to. 



the bolometric equivalent energy and tlierefore neutrino production; c/., Eq. (Pj)|j We have 
checked that this conclusion does not depend on the star formation rate or the correction 
factor E{^z) , even if one of the other star formation rates given in Ref. 46 is usedn Therefore, 



we only use the distribution in the left panel of Fig. |3]in the following. It is also noteworthy 
that at around z = 1 the number of contributing bursts is quite low (compare to histogram), 
which means that strong fluctuations may be expected in low statistics samples. 

We show in Fig. Ill the neutrino flux (left) and flavor ratio (right) for the quasi-diffuse flux. 
The thick solid curves show the integrated flux, the thin solid curves the contributions from 
individual redshift ranges (left panel only). The dashed curves show the reference values 
from Table [TJ Obviously, the diffuse flux peaks at somewhat higher energies than the WB 
reference, because of the main contribution to the flux coming from smaller redshifts. The 
same effect is visible for the flavor ratio. However, neither the characteristic wiggles in the 
shape, nor the flavor ratio transition from pion beam to muon damped source are averaged 
out. The reason is the relatively sharp peak in the weight functions in Fig. |3| which hardly 



^In fact, the IceCube stacking method in Refs. [6|ll| does not compute the isotropic equivalent luminosity, 
needed for the pion production efficiency, back from the observation on a burst-by-burst basis, but uses 
instead a fixed value for the isotropic equivalent luminosity for the pion production. Therefore, mainly the 
break energies of the neutrino spectrum are affected by a different standard value of z. 

^We tested this claim by calculating a correctional factor £{z) for the SFR from Ref. 
the approach from Ref. [45]. Even though this factor leads to more bursts at higher redshifts, the peak 
contribution still comes from bursts at z ~ 1. 
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Figure 5: Neutrino flux (left) and flavor ratio (right) for the reference values from Table fl] where the 
magnetic field B' is varied. The burst with B' = 300 kG is normalized to the WB spectrum, leading to 
normalization factors C' and C' for the input proton and photon spectra. The same proton and photon 
spectra, including normalization factors, are used for the other bursts. 



depends on the star formation rate. One can see from the relative contribution curves 
in Fig. |4] (left panel) that indeed z ^ 1 dominates the flux. The weak dependence on 
the star formation rate is consistent with Ref. 12 . This conclusion obviously changes if 



there are contributions with a different systematics, such as from GRBs originating from 
population III stars 47 , which have different population characteristics. 



3.2 Magnetic field 

Assuming that the magnetic field B' does not significantly change the energy budget of the 
jet, it will only have a significant impact on the maximal proton energy and the cooling 
of the secondaries. We show the effect of different values of B' in Fig. [5J As can be seen 
from the different curves, the flux shape changes due to the magnetic field dependence of 
the synchrotron cooling breaks, see Eq. ^. More precisely, higher magnetic fields lead to 
synchrotron losses dominating at lower energies. Therefore, fewer neutrinos are produced in 
case of high magnetic fields compared to lower magnetic fields. In addition, for B = 7MG, 
the highest peak is no longer determined from the muon break, but the pion break. As 
can be seen from the right plot in Fig. [sj the effect on the flavor ratio R is qualitatively 
different from the effect on the flux shape: depending on the value of B', the transition of 
the flavor ratio from pion beam source to muon damped source shifts in energy, but the 
flavor composition does not change. This shift is also a consequence of the breaks due to 
synchrotron cooling; the higher the magnetic field, the lower the energy of the transition. 
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Figure 6: Distribution of 10 000 bursts as a function of the magnetic field B' (histograms). The left 
panel uses the log- linear distribution with the 3cr range of the exponent going from 70 kG ^ B' < 700 kG, 
the right panel shows the broader version with the range B' = 7kG . . . 7MG, see main text for details. 



Therefore, it is obvious that the flavor ratios may be a direct handle on the magnetic fleld 
inside the source. Note that for B = 7MG, the small dip in the flavor ratio at low energies 
is the result of a pile- up effect. 

For the magnetic fleld, the main problem is the theoretical distribution of the values of B', 
which is unknown. First of all, note that the value of B' enters linearly in Eq. Q, describing 
the spectral break from the energy losses. On the other hand, it depends only on the square 



root of the energy budget, see Eq. (13). Therefore, strong variations of B' will drastically 



affect the spectral shape, whereas the energy budget is hardly affected. In order to get an 
educated guess for the distribution of B', one may vary e^/ee from 0.1 to 10 with a normal 
logarithmic distribution, leading to 70 kG ^ B' < 700 kG in the FB-S (or CB) case. In 
the FB-D case, one can easily see from Eq. (17) that variations of ii^iso,boi, tv, Tgo, and z 



also have only a moderate impact on B', whereas F enters to third power. For example, 
for r = 100 . . . 1000, we have B' = 7kG . . . 7MG, suggesting a broader distribution of B'. 



This, however, depends on the view point: it is correct if the observables in Eq. (17) are 



used as input, and the quantities are calculated back to the SRF. However, if the bursts 



instead share similar characteristics in the SRF, such as E'^^^^^^i, N, and V^g^, see Eq. (13), or 



even -Eiscboi, Re, and Ad' (these parameter sets are related by F, but not B'), the narrower 
distribution is more suggestive. We therefore choose two distributions, which are shown in 
Fig. |6] and deflned as 



B' 

with P(x) 



10^ G 



OC e 2a2 



(20) 



16 



10" 



5 10" 






■^ 10" 



10" 



NeLiCosmA2()ll 




Standard spectram ^^"v^ 




^^^^^^^ 




if 


W D-wide 


If 


uV /^ 


1/ 
1/ 

It 
1/ 
1/ 
1/ 
It 
if 
It 


|W^ D-narrow 


It 


. "X 


It 
It 
It 
1/ 


\ \ 


11 


l\ 


11 


^ \l » 


It 


^ \l 




\ ll 


1 1 \ >\J 



lO'' lO'* lO' 10** 10' 10* lO' 

fy/GeV 



0.8 



0.7 



0.6 



<« 0.5 



0.4 - 



0.3 - 



0.2 



NeuCosmA20ll 

Muon damped 


WB plateau 
Z)-narrow 


- 




V ?i 

If J 
lijr 


D-wide 

Pion beam 




1^ 


"Neutron beam 




- 






" 



£„/GeV 



Figure 7: Neutrino flux (left) and flavor ratio (right) for the magnetic field distributions in Fig. ^ 
The thick solid curves show the flux for the two magnetic field distributions. For comparison, we show the 
reference burst with values from Table [l] as a dashed curve. The thin dashed curve in the left panel is the 
WB reference, the total fluxes are normalized to. 

Here, x is the center of the peak in a;, and a is the standard deviation of the Gaussian. For our 
simulations we used a central value x = 5.34, which corresponds to B' ~ 289 kG (~ 300 kG). 
The standard deviation for the first distribution (narrower distribution, left plot of Fig. IgI) 
is 0"! = 0.17, while the standard deviation for the second (wider distribution, right plot of 
Fig. Q isa2 = 0.500 

We show the diffuse flux from 10 000 individual bursts in Fig. [7} using the magnetic field 
distributions in Fig. [6j In both panels, we can see that the results for the narrow distribution 
(D-narrow) do not deviate much from the results of a single burst (dashed curves). The 
double peak and the contribution from kaons are clearly visible and deviate only slightly 
from the standard burst. However, in case of the broad distribution (D-wide), the flux 
shape (left panel) is washed out to a single peak with the maximum at around 10^ GeV. 
Neither double peak and kaon contribution are visible anymore, which comes from the 
different break energies for the different bursts because they depend on the magnetic field 
B'. Hence, the neutrino flux below the (lowest) muon break is unaffected by the magnetic 



^°Note that the Gaussian distribution of the exponent leads to a slightly different distribution on a linear 
scale. Most notably, the peak is no longer at the median value of the Gaussian distribution, but at a smaller 
value. This comes from the conservation of probabilities during the change from logarithmic to linear scale 
and leads to an additional suppression of high B' values with B'~^ on a linear scale, see Fig. [6] In addition, 
note that we assume that a change of the magnetic field does not effect the amount of energy present in 
protons and photons. However, if the total energy was fixed, a higher magnetic field would also shift the 
energy fractions, leading to less energy available for protons and photons, and therefore even less neutrino 
production. 
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field distribution. For the flavor ratio (right panel), the transition from pion beam to muon 
damped source occurs over a wider energy range for the distribution "D-wide" than for 
"D-narrow". As a consequence, R carries two pieces of information: the (mean) magnetic 
fleld and the broadness of the magnetic fleld distribution, while being hardly affected by 
any other astrophysical parameter, radiative process, etc. (apart from F, which shifts the 
energy of the transition). This information can then be translated into clues for the models, 
such as how alike the bursts are actually in the SRF. 

3.3 Lorentz boost 



As we have discussed in Sec. 2.2, one of the main issues for the Lorentz boost F is the 
translation of F into the relative contributions of the bursts, which is model-dependent. 
Here, we test three different hypotheses: 



i oivr . -Ti-b a uuiiacquciiuc, j^ ^ '^^ • 

ability to observe a burst is oc 1/F^. 



CB Cannonball-like model, bursts alike in SRF. As a consequence, E^cf)^, oc F . The prob- 

^2 



FB-S (Fireball-Shock) Fireball model, bursts alike in SRF. As a consequence, -E^^jy oc 
F^. The probability to observe a burst is oc 6'^ [9: jet opening angle). 

FB-D (Fireball-Detector) Fireball model, bursts alike at the detector. As a conse- 
quence, El(f)y OC F^"^. The probability to observe a burst is oc 6'^. 

These three interpretations may include the most extreme cases, which means that the 
actual consequences of different values of F may lie somewhere in the middle. Note that 
assumptions close to hypothesis FB-D have been typically implied in the literature, such as 



Refs. 17, 18,48], whereas Ref. 21 points towards FB-S 

We show in Fig. Islthe impact of F on the neutrino flux (left) and flavor ratio (right) for the 
GRB with the reference values from Table [T] for hypothesis CB. In the shown plots we have 
flve different values to illustrate the effect of the Lorentz factor. The curve with F = 10^'^ 
represents our standard burst from Table [T| while the values of F = 100, and 1000 are 
commonly used minimal and maximal values - although Fermz-LAT has observed GRBs 
with higher values of F, such as for GRB 080916C. From Eq. ([s]), it is clear that bursts with 
large Lorentz factors dominate, which also shift the flavor ratio transition from pion beam 
to muon damped source to higher energies. The strong Lorentz factor dependence comes 
from the Lorentz boost of the energy and the forward collimation of the Lorentz cone for 
CB: the opening angle is proportional to 1/F^. 

Although the theoretical distribution of F is unknown, there exist empirical analyses. For 
the sake of simplicity, we assume that it should roughly reproduce the observed (derived) 
distribution in Ref. p] (c/.. Fig. lb)FH Therefore, we empirically use a probability distribu- 
tion 
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(logio(r-50)-;E)2 1 
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o j-> 


F-50 



n^)^^.^ — -^ — T^^ ' (21) 



^^Note that this analysis is based on the maximally observed photon energy, whereas different approaches, 
such as based on the afterglow, may be more reliable; see, e.g., Ref. [49]. However, while the peak contri- 
bution in r may be somewhat affected by the distribution function, the qualitative results from this section 
depend rather on the tested hypothesis, i.e., interpretation of different values of T. 
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Figure 8: Neutrino flux (left) and flavor ratio (right) for the reference values from Table fl] where the 
Lorentz boost T is varied. The burst with F — 10^'^ represents our standard burst as in Table [l] and is 
normalized to the WB spectrum. In the shown figures the hypothesis "CB" is used, which determines the 
scaling of the flux with F, i.e., the proton and photon spectra arc normalized to the standard burst, and 
only F is varied among the bursts. 
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Figure 9: Distribution of 10 000 bursts as a function of the Lorentz factor F (histograms) and the 
relative contribution of a value of F to the total diffuse flux (thick solid curves) . The relative contribution 
consists of the distribution of bursts and a scaling factor depending on the hypothesis (see main text). 
Additionally, the dashed thin curves represent the analytic version of the distribution of the number of 
bursts. The different panels correspond to the different hypotheses, as explained in the main text. 
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where the last factor comes from the assumption of a log normal distribution. In the cases 
FB-S and FB-D, the jet opening angle enters the probability distribution, which may be 
correlated with F. We use n = 1 with the median value x = 2.6 and the standard deviation 
a = 0.25 to reproduce the observed distribution in Ref. [s], which may be described, for 
instance, by ^ oc 1/vT. Note that exactly this relationship has been shown to be consistent 



in Ref. 21 with the Ghirlanda relation [50] for the hypothesis FB-S. In the CB case, we use 
instead n = 2 with x = 2.7 and the standard deviation a = 0.25, which can be theoretically 
motivated by the probability oc 1/F^ to observe a burst, coming from the forward coUimation 
of the Lorentz cone 



im 



In Fig. [9] we show the distributions for all of the three hypotheses. In each plot, the 
histogram shows the actual number of bursts generated with our Monte Carlo generator, 
while the thin dashed line is the theoretical distribution of the number of bursts depending 



on Lorentz factor as in Eq. (21). The thick solid curves represent the relative contributions 
to the diffuse flux, which are the respective distribution functions times the appropriate 
scale factors for the models, as described above. Even though all three histograms look 
similar by construction, the relative contributions exhibit significant differences. While the 
main contribution for the CB model comes from bursts with F ^ 700 — 750, the fireball 
models have their main contribution from lower values. The main contribution in the FB-S 
model comes from bursts with F ^ 400, which is still above the peak of the F distribution. 
For the FB-D model it is even lower with F ^ 140 below the peak of the F distribution. 
In none of the three cases the peak of the number distribution corresponds the peak in the 
relative contribution. Note that we have also checked that an intrinsic distribution of F, as 
it may be already expected within one object, such as a relativistically expanding fireball 
from which different components are observed with different viewing angles and therefore 
Doppler factors, is narrower than the distributions used here. In addition, note that the 
contribution functions of CB and FB-S are qualitatively similar. For CB, the probability 
to observe a burst scales oc 1/F^ and the flux £"^0,^ oc F"^, i.e., a total factor of F^ survives 
in the contribution function. For FB-S, if 6 is fixed, one has a flux El(f)^ oc F^ with the 



same scaling^ As a consequence, which can also be seen in Fig. |9| for CB only few bursts 
contribute at the peak of the contribution function, whereas in the case FB-S many bursts 
contribute. Therefore, generically, one expects large fluctuations in stacked neutrino fluxes 
in the CB case, whereas in the fireball cases the fluctuations from F are expected to be 
smaller. We will address this issue in Sec. HI 



The resulting diffuse neutrino fluxes for the three hypotheses can be found in Fig. 10 One 



can see from the left and the middle panels of Fig. 10 that the resulting fluxes and flavor 



ratios for CB (left panels) and FB-S (middle panels) are quite similar. Both show the 



distinct features of a single burst we already discussed in Ref. 14 . The total fluxes (thick 
solid curves in upper panels) are just slightly shifted in energy because of the different 
scaling with Lorentz factor F. The effect of the different contribution functions is depicted 
by the thin solid curves representing the contributed flux from all bursts with F < 300, 
< 500, and < 1000, as labeled in the plots. The contribution from high F bursts in the CB 



^^Note that the offset of 50 is included to reduce the number of bursts with F < 100. 
^■^In order to reproduce the empirical distribution, we effectively use 9 oc l/vT, which explains the 
difference in the contribution functions. 
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Figure 10: The diffuse ffuxes of muon neutrinos (upper panels) and the flavor ratios (lower panels) from 
10 000 accumulated bursts. The upper plots show the contributions from all bursts with Lorentz factors 
r < 300, < 500, and < 1000 separately. The "All T" curves are normalized to the WB spectrum. The 
panels from left to right show the CB model, the FB-S model, and the FB-D model. For comparison, we 
show the reference burst with values from Table [T] as a dashed curve. 



model is significantly higher than in the FB-S model because the peak contribution for CB 
comes from burst with F ^ 700 — 750, while for FB-S, it comes from bursts with F ^ 400. 

For FB-D (right panels), the results are qualitatively very different. This comes from the 
assumption of fixing all parameters at the detector, not in the SRF, as in the other two 
models. Thus, for each burst we have to calculate the isotropic volume with the correspond- 



ing collision radius and the corresponding shell thickness, see Eqs. (15) and (16). Moreover 



all energies have to be boosted to the SRF with the correct Lorentz factor. In all burst 
simulations in the FB-D case, we scaled C' C oc F" 



^7,break OC F \ and B' (XV 3, c/., 
Eq. (17). In consequence, the variation of the normalization factors leads to a dominance 
of low F bursts with F ^ 140. Additionally, the change of the photon break energy in the 
SRF e^ break ^'^d the magnetic field B' lead to a change of the flux shape. Thus, the features 
of a single burst are smeared out by summing over different bursts, and the flavor ratio 
transition broadens significantly. Note, however, that the smearing of magnetic field effects 
is the result of an implicitly assumed wide distribution of B' with a similar effect as in the 



previous section (which was also estimated from the F dependence of Eq. (17), applicable 
to FB-D). 
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Figure 1 1 : Neutrino flux for the reference values from Table n] where the lower spectral index a-,, (left 
plot), the upper spectral index /3^ (central plot), or the photon break energy e' break (right plot) of the 



target photon spectrum is varied. The spectra labeled a^ 



-1.00, /3t, 



-2.00, and e: 



7, break 



= 14.76 keV, 



respectively, represent our standard burst as in Table [T] and are normalized to the WB spectrum. The 
other bursts are normalized by assuming that the total energy in photons is the same. 



In summary, the effects of a F distribution are model-dependent. Wliile most of tlie literature 
focuses on case FB-D, where magnetic field effects of the secondaries are smeared out, 
different results may be expected even within the fireball model if the bursts are similar in 
the SRF. Since the observational data may not be mature enough to discriminate among 
the different cases, a distribution of F has to be interpreted with care, and cannot be 
performed at a similar level as a distribution of redshift. On the other hand, it means that 
some of the high-F bursts accessible by Fermz-LAT could indeed contribute more to the 
neutrino flux than previously anticipated. The test of magnetic field effects, such as the 
multi-peak structure of the spectral shape or a sharp transition of the flavor ratio, may 
help to discriminate FB-D versus CB and FB-S. In particular, the observation of these 
features may exclude conventional fireball phenomenology (FB-D). However, CB and FB-S 
are intrinsically hardly distinguishable. Perhaps the statistical fluctuations of low statistics 
samples may give some clues, see discussion in Sec. |4} 

3.4 Photon spectrum 



In the following, we test the effect of the photon spectrum on the neutrino flux, i. e., vary the 
parameters in Eq. (|2|). In Fig. 11, the effect of a.y (left plot), /3^ (middle plot), and £lj,_break 
(right plot) on the v^^ flux shape are shown. As before, the spectra with the standard 
parameters taken from Table [T] are normalized by matching the integrated muon neutrino 
energy to the one in the WB flux. For the subsequent calculations we, however, fix the 
total energy in photons by integration over the whole spectrum. Before we come to the 
effect of the individual parameters, note that increasing the spectral indices or decreasing 
the break energy makes the spectra softer because of this normalization constraint. This 
means that more target photons for high energy protons are available, and the photohadronic 
interactions and therefore the neutrino production becomes enhanced. Although the photon 
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Figure 12: The distributions of the spectral parameters, namely the first spectral index a^ (left 
plot), the second spectral index /3^ (middle plot), and the photon break energy e' ^i-oak (right plot). The 
histograms depict the 10 000 randomly picked values following the distributions obtained in Ref. fsTl. The 
dashed curves represent the analytic version of the distribution function. 



spectrum may, in practice, be correlated with the other burst parameters in a non-trivial 
manner, this approach will turn out to be sufficient to demonstrate our main points. 

In the left and middle panels of Fig. 11, we show the effect of a^ and /3^, respectively. As 
expected, the softer the spectrum is, the more neutrinos will be produced. For a^, the 
neutrino spectral index above the ffist break varies, for /3^, the neutrino spectral index 
below the ffist break, where these variations are anti-correlated with the photon spectral 
indices, as explained in Sec. |A.2| For el, break (I'ight panel), a lower break energy in the 



'7, break 

photon spectrum leads to a higher -E,y,break- The effect is, however, hardly visible on the 
left end of the spectra, because the muon and pion cooling breaks dominate the spectral 
shape, which are (in all cases) not shifted in energy by the photon spectrum. Only the 
relative contributions between muon peak and pion peak may change with e^ break; because 
the relative importance of the multi-pion processes is altered. 

We now want to introduce the distribution of the different parameters. Compared to the 
often used standard values a^ = —1 and /3^ = —2, we use data of long duration GRB 
obtained in Ref. |5l], based on Fermz-GBM data, to model our distributions, which are 
shown in Fig. 12 The distribution for a^ peaks at a^ = —0.90, is Gaussian, and has a 



standard deviation of a^ = 0.36. Therefore, our standard value is still in the la range 
of the measured peak from that catalogue. For /3^, there are only 84 long bursts in the 
data catalogue with a time-integrated spectrum best fit with the Band function. We have 
fitted a Gaussian to these values with a mean value /3^ = —2.35 and ap = 0.16. Here, our 
standard value is only within the 3a range of the distribution. However, the statistics are 
considerably lower compared to the amount of data used to obtain the distributions for a^ 
and £l),_break- For e^break ^^ kcV, the peak of the distribution is on a logarithmic scale at 0.18 
(lO^'^^keV = 1.51 keV) with a standard deviation a = 0.36. As one can see the observed 



break energies are significantly lower than our standard value of el 
obtained by assuming -Ej/^break 



7, break 



14.76keV we 



10 GeV. Thus, the observed bursts should have higher 
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Figure 13: Neutrino flux (left plot) and flavor ratio (right plot) for a diffuse muon neutrino flux obtained 
by summing over 10 000 individual bursts with varied parameters of the input photon spectra. For the thick 
solid curve, denoted as "All varied", both photon indices, a^ and /3^, as well as the break energy e' 



7, break 



are 



varied in every spectrum. However, the total energy in photons is considered constant and each spectrum 
is normalized accordingly (by integration). In the left plot, the thin solid curves represent the resulting 
(quasi-) diffuse flux if only a^, (3^, or e^ is varied. The curves are labeled accordingly. The thin dashed 
curve represents the WB spectrum, the dashed curve the standard spectrum we obtained earlier. 



values of ^i.,breakQ 

The resulting neutrino fluxes and flavor ratios from the distributions in Fig. 12 are shown 
in Fig. 



13 



in the left and right panels, respectively. Here, the parameters a^, /3.y, and £^^ break 
have been varied separately (thin dashed curves) and simultaneously (thick solid curves). 
As the main result, the flavor ratio is hardly affected by the photon spectrum variation, and 
the generic structure of the flux shape also does not change. The most signiflcant change 
may come from the different mean value of £^ break compared to our reference burst, and 
the relative enhancement of the second peak coming from pion decays. These parameters, 
however, do not have any impact on our line of argumentation. The reason is that the peaks 
in the SRF are determined by Eq. (Il]), which is completely independent of the shape of the 
target photon spectrum. 



14 



Note that in Ref. 



51 



the distribution for the observed Epcak is given, and not for the break energy in 



the SRF, e' break- Therefore, we must boost the distribution from Ref. 51 to the shock rest frame under 
consideration of our standard values for the Lorentz factor F and the redshift z, see Table [ij Here, we 
neglect the slight difference between -Epeak from the Band function parameterization 52 and e^, break of the 
broken powerlaw parameterization in Eq. (|2|). 
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Figure 14: Example for a smaller subset of bursts. In the left panel, one sample of 100 randomly picked 
bursts for the strong evolution case is shown (histogram), for which the distribution function is illustrated 
as dashed curve. The solid curve is the (theoretical) contribution of a redshift value to the diffuse flux (scale 
arbitrary). In the right panel, the thick solid curve shows the resulting quasi-diffuse flux for this sample. 
This has to be compared to the thick dashed curve, which refers to the high statistics result for the variation 
in redshift, see Fig. l4J The other (thin) solid curves show the result for other subsamples, where the curves 
labeled "high sample" and "low sample" show extreme possibilities obtained in about 0.1% of all cases only. 
Moreover, the classic WB flux shape is shown as a thin black dashed curve. 

4 Quasi-diffuse fluxes from low statistics samples 



Here, we discuss the aggregation over a finite number of bursts, as it is used in a stacking 
analysis, see, e.g., Ref. [6]. We use tlie same metliod as before, but we vary the discrete 
number of bursts n used to compute a quasi-diffuse flux. Obviously, the statistical fluctua- 
tions coming from the summation over significantly fewer bursts than in a diffuse flux will 
introduce a systematical error in the extrapolation from the burst sample to a quasi-diffuse 
flux. However, it is exactly this extrapolation which is needed to compare the limits from 
neutrino telescopes to theoretical predictions, such as the Waxman-Bahcall estimate 30 



and to compare to predictions from fireball phenomenology. In this section, we quantify 
this systematical error using the redshift distribution, since we have shown that conclu- 
sions based on redshift may be drawn in the least model-dependent way. However, we also 
demonstrate how statistics can be used to discriminate among different hypotheses for F. 



In Fig. 14 we show one example for a sample of 100 randomly picked bursts (left panel). From 
the comparison between the histogram and the contribution function (solid curve), which 
depicts the weighted contribution of individual bursts, it is clear that strong fluctuations can 
be expected. The resulting quasi-diffuse flux for this distribution is shown in the right panel 
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Figure 15: Probability that the quasi-diffuse flux extrapolated from a low statistics sample with n 
bursts is larger than a certain fraction of the diffuse limit (see legend for different values of n). This 
function corresponds to (one minus) the cumulative distribution function of the probability density. Here, 
the strong evolution case for the redshift is used for the burst distribution function. 



(thick solid curve), wliicli somewhat deviates from the diffuse high statistics limit (dashed 
curve) because of these statistical effects. In addition, two extreme cases of samples of 
100 bursts are shown as thin solid curves; such extreme cases are only obtained in about 
0.1% of all cases. If other (independent) parameter distributions are present, such as for 
B' or r, we expect that the resulting statistical error, to a first approximation, adds up 
in a Gaussian manner. Therefore, the redshift distribution will only give a lower estimate 
for the introduced systematical error. A more detailed discussion including the impact of 
the gamma-ray detection threshold on a stacking analysis using a simultaneous redshift- 
luminosity variation can be found in App. [B| 

While this figure illustrates that some effect of the statistical fluctuations can be expected, 
which depend on the sample size, it is not sufficient for a quantification of the problem. 



Therefore, we choose the following method: We show in Fig. 15 the probability that the 
resulting quasi-diffuse flux from a low statistics sample with n bursts is larger than a certain 
fraction / of the diffuse limit, ^^ By definition, the probability that the diffuse flux is larger 



than the fraction / < 1.0 of the diffuse limit is always one, whereas it is zero for / > 1.0, 
which is shown as the step function labeled "diff. limit" in the figure. The smaller the 
sample size is, the stronger is the deviation from this curve. For example, for tt, = 5, the 



^'^Methodologically, we first systematically generate bursts with < z < 6 in steps Az = 0.001 and 
calculate the energy flux density. Then, in the second step, we draw samples of n bursts with probabilities 



according to Eq. (19 1 (with strong evolution). In order to compute the probabilities in the figure, we 
repeat this procedure for each sample size N times, where, ideally, N —^ oo, and, practically, large enough 
N = 100 000 are used. The used energy flux density ^ E^ 4>ti{Eu) d-Ejy is obtained by integration performed 
from 10° GeV to 10^° GeV. 
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n 


±10% 


±20% 


±50% 


5 


0.07 


0.14 


0.40 


10 


0.09 


0.18 


0.52 


50 


0.14 


0.30 


0.82 


100 


0.17 


0.37 


0.91 


300 


0.23 


0.50 


0.96 


1000 


0.30 


0.69 


0.98 


10000 


0.48 


0.97 


0.99 



n 


Rel. error 90% CL 


Rel. error 3cr 


5 


0.15-2.45 


0.09-20.89 


10 


0.23-2.22 


0.13-19.05 


50 


0.44 - 1.72 


0.30 - 10.28 


100 


0.53-1.57 


0.39-8.78 


300 


0.64-1.38 


0.53-6.53 


1000 


0.72-1.25 


0.64-5.15 


10000 


0.83-1.08 


0.78 - 2.62 



Table 2: Lower row, left table: Probability that the quasi-diffuse flux is within a certain interval of 
the diffuse limit for different stacking sample sizes n. Right table: Systematical error coming from the 
extrapolation of the quasi-diffuse flux from a sample with n bursts at different confidence levels. The 



corresponding figures (upper row) illustrate how this information is obtained from Fig. 14 for n = 100. 
Here, the strong evolution case for the redshift is used for the burst distribution function. 



probability that / < 0.2, i.e., at most 20% of the diffuse limit are obtained as a result, 
is about 1 — 0.89 = 0.11, and the probability that / > 5.0, i.e., five times more than the 
diffuse limit is obtained, is about 0.01. While this figure is somewhat difficult to read, it 
contains already all of the information discussed in the following. Note that this method 
corresponds to computing one minus the cumulative distribution function of the probability 
density. 

For example, one can ask what the probability for a quasi-diffuse flux, derived from the 
stacking of n bursts, is to be within a certain fraction of the diffuse (high statistics) limit. 
This question is addressed in the left column of Table [2j The plot in the upper left panel 
illustrates how this information is obtained from Fig. [l4j for a certain interval around the 
diffuse limit curve, the corresponding probability contained in this interval is read off from 
the vertical axis. From the table, one finds that n = 10 000 is already close to the diffuse 
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limit. Quantitatively, the probability that a quasi-diffuse ffux based on 10 000 bursts, such 
as based on the number of observable bursts during ten years of IceCube operation, is 
within 20% of the diffuse limit from a purely statistical point of view, is 97%. That is 
sufficiently diffuse for any practical purpose, but, nevertheless, it is surprising that small 
fluctuations are expected. Practically, it means that the diffuse flux during ten years of 
IceCube operation may not be exactly the same as the one during the following ten years at 
the level of 20%. One can also read off from that table that for n = 100, as used in recent 
IceCube stacking limits, the probability to be within 20% of the diffuse hmit is only 37%. 
Only if n = (9(1000), reliable predictions from stacking analyses can be expected. 

A different but complementary question is addressed in the right column of Table |2] Here, 
the variation of the flux with a certain probability window is discussed, see upper right 
figure for an illustration. Since the curve shown in this figure corresponds to the integrated 
probability density, it can be used to derive the range (horizontal axis) in which 90% of 
all cases and 99.73% of all cases are found. These numbers can be used to estimate the 
systematical error in the diffuse-flux extrapolation at the corresponding confidence level. 
The result is shown in the lower right panel in Table |2j It can be inferred from this table 
that the systematical error in the calculation of the quasi-diffuse flux in terms of the energy 
flux density (and thus, at least roughly, in terms of the flux normalization) is at the level of 
about 50% for 100 bursts, 35% for 300 bursts, and 25% for 1000 bursts (90% CL) from the 
redshift distribution only - additional parameter variations will lead to a larger effect (see 
App. [B]). This is to be compared with the naive statistical estimate if all bursts were at the 
same z: the la (90% CL) relative error for 100 bursts may be estimated as 1/v^lOO ~ 0.1 
(0.16), to be compared to the actually obtained 50% at the 90% CL. 

Apart from the discussion of finite sample sizes, the special case n = 1 is interesting. In 
particular, what is the probability that a single bright burst is detected in ten years of 
IceCube operation, given the redshift distribution (strong evolution case) only? In this 
case, the multi-messenger connection to gamma-rays may not be necessary, and the burst 
should also outshine the backgrounds. As a prerequisite, a few events from such a single 
burst, say at least three, have to be detected. The probability for such a burst depends on 
the level of the diffuse flux. If it saturates the WB estimate, such a burst would need to 



contribute only about 1/30 of the total flux If the diffuse flux is only at the level of 1/10 of 



the WB estimate, such a burst has to contribute 1/3 of the total flux. From these numbers, 
we can compute the maximally allowed redshift, which is -Zmax.i — 0.14 and 2;max,2 — 0.05, 
respectively, for these two different cases. As the next step, we obtain as probabilities for a 
single burst P{z < 0.14) = 5 • 10"^ and P{z < 0.05) = 2 • 10"^ Given that C(10 000) bursts 
are expected in the observable universe in ten years, the probability that such a bright 
burst, which does not require the observation of the gamma-ray counterpart, will occur at 
least once is 1 - [1 - P{z < Zmax)]" with n = 10 000, which is about 40% if the WB flux is 
saturated, and about 2% if the diffuse neutrino flux is an order of magnitude smaller. 



In Sec. |3.3[ we have discussed that the consequences of variations of T are model-dependent. 
However, we have illustrated that the models CB and FB-S lead to relatively similar fluxes 
in the quasi-diffuse limit, whereas for low statistics, stronger fluctuations are generically 



^^This can be derived from the estimate that the WB neutrino flux, if saturated, leads to 0(100) events 
in ten years of operation of IC-86, which is consistent with the current IC-40 bound. 
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Figure 16: Figure similar to Fig. 15 but for a variation of the Lorentz factor F for the cases CB and 



FB-S. Here, the two models are compared for n = 50 and n = 300. 



to be expected for CB than for FB-S (c/., Fig. |9| wliere tlie bursts statistics close to the 
peaks of the distribution functions are different). We test this hypothesis quantitatively in 



Fig. [16} where we directly compare the integrated probability distributions between CB and 
FB-S for two sample sizes (n = 50 and n = 300). It can be clearly seen that the case CB 
deviates much stronger from the diffuse limit (step function), as expected. In addition, this 
difference between the two models becomes smaller for larger n, i.e., in the high statistics 
limit, there is hardly any difference. For example, for n = 50, the probability to be within 
±10% of the diffuse limit is only 10% for CB, whereas it is about 40% for FB-S. This 
statistical effect suggests to use this effect for a test of the model: For example, one may 
work with randomly chosen n = 50 samples from a larger stacking sample accumulated over 
many years. Then by the test of the variance of the neutrino fluxes among these samples 
one may draw conclusions on the underlying model. Note that CB/FB-S versus FB-D can 
be discriminated by the impact on the magnetic field effects, as we have discussed earlier. 



5 Summary and conclusions 

We have numerically computed ORB neutrino fluxes and flavor ratios including different 
meson production modes, magnetic field effects, and flavor mixing. Our starting point 
has been the minimal set of assumptions needed to compute the neutrino fluxes at the 
same level as it is used in recent IceCube stacked limits, but including the full spectral 
dependencies and the cooling of the secondaries explicitely. We have kept these assumptions 
as model-independent as possible. As the next step, we have connected the quantities 
needed for the neutrino production to different models, such as a cannonball-like model and 
conventional fireball phenomenology, where we have clearly demonstrated where model- 
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dependent assumptions or relationships enter. 

The main purpose of this study has been the computation of aggregated fluxes, where we 
have distinguished diffuse fluxes and quasi-diffuse fluxes obtained from stacking analyses. 
For the computation of the aggregated fluxes, we have computed the flux burst by burst, 
where we have used 10 000 bursts for the diffuse flux - as may be expected in the observable 
universe within ten years of neutrino telescope operation. We have established that these 
10 000 bursts lead to a flux close enough to the diffuse limit, i.e., that the systematical error 
coming from the burst statistics is small in that case. Quantitatively, from the redshift 
distribution (strong evolution case), the neutrino flux in these ten years will be within 20% 
of the diffuse limit with 97% probability. 

Our computations of aggregated fluxes have been based on parameter distributions for 
redshift z, magnetic field 5', Lorentz boost F, and the spectral shape of the target photons. 
We have identified two major issues: The theoretical distribution of the parameter, and the 
interpretation of a different parameter value in terms of the relative contribution to the flux. 
A relatively simple case is redshift, which may follow the star formation rate in some form, 
and the interpretation of the relative contributions of the bursts, if assumed to be otherwise 
comparable at the source, is straightforward. A completely different case is F for which 
the theoretical distribution is unknown and the translation of the relative contribution of 
the bursts is model-dependent. For example, in the conventional fireball model, the size of 
the interaction region is estimated from the observed variability timescale t^, which implies 
that it strongly depends on the Lorentz factor F. As a consequence, small values of F 
are assumed to dominate the neutrino flux. In the same model, the hypothesis of similar 
properties in the comoving frame leads to the dominance of high F, as it is consistent with 



the recent Ghirlanda et al. study 21 . As a consequence, some of the Fermz-LAT accessible 
bursts with high F may contribute stronger to the neutrino flux than previously anticipated. 
Since the underlying assumptions for the different parameter distributions are qualitatively 
different, we have not mixed these different cases in the computation of the diffuse fluxes; 
in addition, our quantitative statements (see below) are based on the redshift distribution 
only, as the most conservative assumption. 

The synchrotron cooling versus decay of the secondaries (pions, muons, kaons) in combina- 
tion with multi-pion photohadronic processes leads to a characteristic multi-peak structure 
of the GRB neutrino flux. In addition, a transition of the flavor composition from a pion 
beam source (neutrinos from pion decay chain) to a muon damped source (neutrinos from 
muon decays suppressed) at high energies is expected, which can, in principle, be measured 
by the muon track to cascade ratio in a neutrino telescope. We have established that these 
features are, in general, also present in diffuse fluxes, such as if the redshift of the bursts 
is varied. Only too strong variations of the magnetic field may destroy the spectral peaks 
in the flux and smear out the flavor ratio transition. In this case, the spectral shape of 
the flux cannot be used anymore to obtain information on the magnetic field, whereas the 
flavor composition depends on the mean value of B' and the broadness of its distribution 
with no other significant parameter dependencies. This means that the flavor ratio can be 
regarded as the most robust prediction one can make for diffuse GRB neutrino fluxes. The 
observation of magnetic field effects can also be useful as a model discriminator: we have 
demonstrated that the observation of spectral features in diffuse GRB neutrino fluxes or 
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a narrow transition of the flavor ratio can discriminate between the conventional flreball 
assumption that bursts with low F dominate the neutrino flux and the hypothesis that the 
bursts are alike in the comoving frame. 

Apart from diffuse fluxes, we have also discussed quasi-diffuse limits obtained from the 
stacking of a relatively small number of bursts, such as C(IOO) bursts used in recent IceCube 
analyses. The stacking statistics of a small number of bursts has been shown to introduce a 
systematical error in the computation of the quasi-diffuse flux, which is at the level of about 
50% for 100 bursts, 35% for 300 bursts, and 25% for 1000 bursts (90% CL) from the redshift 
distribution only - additional parameter variations will lead to a larger effect (see App. [B| . 
Therefore, (9(1000) bursts may have to be stacked for reliable conclusions. The reason is 
that the redshift z leading to the largest contribution to the flux does not coincide with the 
maximum of the redshift distribution. As another consequence, in stacking analyses, it may 
also be consistent to assign z ~ 1 to a redshift not measured for a burst, since, depending on 
the method, a different value may lead to an overestimate of the neutrino flux or to wrong 
energies of the spectral neutrino breaks. This conclusion holds almost independently of the 
redshift distribution, including the strong evolution case. Similar standard values can be 
inferred for F, which are, however, model dependent. For example, F ~ 700 gives the main 
contribution in our cannonball ansatz, F ~ 400 in the flreball model with the parameters 
flxed in the comoving frame, and F ~ 140 for the conventional flreball assumption. While 
normally a high stacking statistics is desirable, we have also demonstrated that indeed 
the fluctuations among different stacking samples may be used as a discriminator between 
cannonball and flreball approaches. 

A very interesting special case is the observation of a single bright burst detectable without 
multi-messenger counterpart. From the redshift distribution (strong evolution case), we 
have estimated that the probability to observe at least one single bright burst individually 
detectable in IceCube during ten years of full-scale operation is about 40% if the Waxman- 
Bahcall flux is saturated, and 2% if the diffuse neutrino flux is an order of magnitude smaller. 
Thus, at best, there may be a 50-50 chance for such a single burst detection. 

We conclude that neutrino observations of GRBs, both in terms of flux and flavor ratio, 
may not only help to test the paradigm if GRBs are the sources of the highest energetic 
cosmic rays, but also give important clues on the underlying models. Especially spectral 
features and the flavor composition turn out to be very interesting discriminators of very 
basic assumptions if a signal is observed. If no signal is observed, a systematical error coming 
from the stacking statistics has to be regarded in the context of the limit calculation. 
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A Details of the source model 

In this appendix, we describe details of the source modeL 

A.l Components of the model 

In the following, primed quantities refer to the shock/bulk rest frame (SRF), and unprimed 
quantities to the observer's frame. In some cases, which we denote explicitely, the unprimed 
quantities may also refer to the source (engine) frame. 

Photon and proton input spectra 

We assume that the target proton field in the SRF is given by Eq. ^. Here, the maxi- 
mal proton energy E' ^^^^ is derived from the comparison between acceleration rate of the 
protons 53 

, , c^ qB' 

Cc = ^-^ , (22) 



and the synchrotron loss rate 
leading to 



e^ 5'2 E' 
9tc Eq m'^c' 



t'syliE') = Tr—-ir. ^ (23) 



1 



E;,_^2-10«(^)^ I,,^) GeV , (24) 



!LY ( B' y 

0.1/ VlO^G^ 

where rj is an acceleration efficiency. This assumption is somewhat "vanilla" in the presence 
of adiabatic cooling (or the dominance of photohadronic cooling), since the maximal proton 
energy may be limited otherwise. However, since the neutrino spectrum will be dominated 
by the energy losses of the pions for E^, > 10^ GeV, we expect these effects to be small. 

The target photon field N'{e') is assumed to be derived from observed GRB photon spectra 
and is parameterized by Eq. ^. The spectral indices a^ and (3y can be obtained from 
GRB observations. The break energy el. break is normally retrieved from the photon spec- 
trum of a burst, but we will need it as a parameter in our model. Moreover, the energy 
^7,max is the maximum photon energy, and e'^^^^^ is the lower energy cut-off of the photon 
spectrum. In practical applications, these may be chosen according to the energy window 
of the gamma-ray observation. Here, we choose e^^nin = 0.2 eV, el^max — ^OOe'^^reak 13 , 
-^p.min — 1 GeV, and ?7 = 0.1. These values are somewhat optimistic in the sense that the 
threshold for photohadronic interactions is always sufficiently passed in the whole energy 
range. However, we have tested that the neutrino results are not very sensitive to e^min 
(as long as e'^^^^^ < l/SOe^break). to e'^^^^^ (as long as e'^^^.^^ > SOe^break). to E'^^^^^ (as 
long as -Ep,min ^ 1000 GeV), and to r] (as long as rj > 0.01). Note that any violation of the 
conditions in brackets for the photon spectrum would be visible even within the BATSE 
energy window (given the Lorentz boost of these quantities). 
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Photohadronic interactions and weak decays 



The accelerated protons interact with the target photon field, where pions, kaons, and 
neutrons are produced. We consider the following processes: 



p + 7 -)■ vr +p , 

p + 7 ^ i\:+ + A/E . 



(25) 
(26) 



Here, p' is either a proton or neutron, vr are one or more neutral, positive, or negative 
charged pions, and A and S are resonances. Note that the effect of kaon decays is usually 
small, but kaon decays may have interesting consequences for the neutrino flavor ratios at 



very high energies, in particular, if strong magnetic flelds are present 54-57 . Therefore 



we consider the leading mode: K^ production (for protons in the initial state) with the 
decay channel leading to highest energy neutrinos. Note that at even higher energies, other 
processes, such as charmed meson production, may contribute as well. 

We follow the description of the photohadronic interactions in Ref. [33] (Sim-B), based 
on the physics of SOPHIA [34]. The treatment, taking into account resonant (including 
higher resonances), direct, and high energetic multi-pion production, allows for an efficient 
computation which is very accurate for power law-like spectra. It predicts well neutrino to 
antineutrino ratios and fiavor ratios. The secondary meson injection rate Q'^^{E'f^ (in units 
of GeV~^ cm~^ s^^) for the interaction between the proton spectrum N'JE') in Eq. O) and 
the target photon field N'{e') in Eq. ([2|) is given by 



Q'.iK) 



dE' , , 



deN'Je')R,{x,y) 



K 



IE' 



(27) 



with X = E^/Ep the fraction of energy going into the secondary, y = {EpS')/mp (directly 
related to the center of mass energy), the "response function" Rb{x,y), and eth the thresh- 
old for the photohadronic interactions (photon energy in proton rest frame). Since many 
interaction types are considered, the response function can be quite complicated; for details, 
see Ref. |3 



From Eq. (27), one can make already two very important observations: First 



of all, in Q'f^{El) only the product of the normalizations of the proton and photon fiuxes 
enters. This means that the ratio between the energy in photons (or electrons) and protons, 
i.e., the baryonic loading, is not explicitely required for the overall normalization of the 
neutrino spectrum - unless processes explicitely depending on the proton or photon den- 
sity contribute significantly [e.g., for inverse Compton scattering). The second observation 
from Eq. (27) is that the steady densities [GeV~^ cm~^] of proton and photon fields in the 
interaction region determine the neutrino production, which we have used as initial input. 
These steady distributions are often {e.g., in Ref. [s]) directly computed from the ejected 
(observable) gamma-ray fiux under the assumption Q'^ ~ A^^/tg^^ (see Eq. (35) below) with 
tggj, related to the size of the interaction region {e.g., shell width in the internaTshock model). 
For the neutrino production, however, it does not matter if the photons escape the source or 
not, which is taken into account in some models such as the photospheric emission model, 
see, e.g., Ref. 58 . For the model-independent part of our calculation, it does not matter. 
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Figure 17: Synchrotron cooling and decay rates for pions, muons, and kaons as a function of energy 
for a magnetic field of B' ~ 300 kG. The gray-shaded region shows the estimated range for the escape 
(or adiabatic cooling) rate if the size of the acceleration region is given by a shell width determined by a 



variability timescale ranging from 0.01s to Is (Lorentz factor F = 300, redshift z — 2; cf., Eq. (15)) 



We let the produced particle species, namely pions, kaons, and neutrons, decay weakly into 
neutrinos as 






v^ + V, 



II ^ 



/i^ -> e~ + z/e + Z//, , 
n — > p + e~ + Ue ■ 



(28) 

(29) 
(30) 
(31) 



The weak decays are treated as described in Ref. 13 , for details see also Ref. 27 . Since 



the decays of muons are helicity dependent, we distinguish between left and right handed 
muons, which has some impact on the flavor ratio. 



Energy losses and escape 

The kinetic equation for the particle spectrum, assuming continuous energy losses, is given 
by (see, e.g., Ref. [59]): 



d ,.,,^,, ,,,,^,,, , N'{E' 



Q\E') = —{h\E')N\E')) + ^^ 
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(32) 



with tggp(£") the characteristic escape time, b'{E') 



-^'OithCUi^') 



■l/E'dE'/dt' 



the rate characterizing energy losses, Q'{E') the particle injection rate (in units of GeV~^ cm"^ s~^) 
and N'{E') the steady particle spectrum (in units of GeV~^ cm~^). Here, we do not consider 
an explicitely time-dependent version of this kinetic equation, because the low statistics in 
neutrino observations will not allow for the resolution of time- dependent features. For the 
charged secondaries (pions, muons, kaons) we assume that they lose energy mostly by syn- 
chrotron radiation, governed by Eq. (|23|, and that they escape mostly by decay 



decay 



mc 



E'To 



(33) 



where Tq is the rest frame lifetime. We show the synchrotron cooling and decay rates for 
pions, muons, and kaons as a function of energy in Fig. [T7| where it can be clearly read 
off that the synchrotron and decay rates are equal at different energies which depend on 
the particle species. The gray-shaded region shows the estimated range for the escape 
(or adiabatic cooling) rate if the size of the acceleration region is given by a shell width 
determined by a variability timescale ranging from 0.01s to 1 s (Lorentz factor F = 300, 
redshift z = 2; cf., Eq. (15)). One can see that synchrotron cooling and decay typically 
dominate, which is an assumption implied in Refs. 5v& . As a consequence, the (steady state) 
spectrum is loss-steepened by two powers above the energy E'^ oc \/rrv^ Jt^ given by Eq. (4), 
where synchrotron cooling and decay rates are equal. Thus, the different secondaries, which 
have different masses m and rest frame lifetimes tq, will exhibit different break energies 
which solely depend on particle physics properties (and the value of B'\ For example, the 
pion spectrum will have a break at 



EL = 1-1 -10' 



5' 



105 G 



-1 



GeV 



(34) 



Note that in the absence of energy losses h{E) = . Then we can easily solve the differential 



equation in Eq. (32) in order to obtain 

N'iE') = Q'iE')t',, 



(35) 



A. 2 Typical results 



Here, we discuss the results qualitatively and quantitatively using the often used Waxman- 
Bahcall flux 



30 as an example; see also Ref. 14 



Qualitative shape of the neutrino spectra 



If the protons are injected with a power law with injection index two, one obtains for the 
prompt GRB neutrino flux, referred to as "WB flux" , 



Ef. 



oc < 



E^ 



-'-'I'.brcak 

E^ 

-'-'f .break 



E^A 



Pu 



E^ 

Ev,w 



for Ey < -Ey,break 

for Eyhj-cak "^ Ey < Ei^T, 



for E^ > E^,^ 



(36) 
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with O;^ = «p — /3^ + 1 ~ +1 and /3j, = a^ — a^ + 1 ~ 0, where this form apphes to the muon 
neutrinos from pion decays. The values for the break energies used in Ref. l30] are -E^,break — 
10^ GeV and E^^.,^ ~ 10^ GeV. For the comparison to the quantitative computation, it is 
useful to define a reference parameter set which reproduces this shape. We use the analytical 
estimates from Ref. [8], assuming that T = 10^'^ and z = 2. The first break energy -E,y,break 
can be related to £-),,break from the threshold of the photohadronic interactions at the source. 
As a minor difference to Ref. [s] , where heads-on collisions between photons and protons are 
assumed for the threshold, we include the effect that the pion production efficiency peaks at 
higher center-of-mass energies (see Fig. 4 in Ref. [33]) to match the numerical results. This 
leads to a factor of two higher photon energy break in the SRF (14.8 keV). The second break 
comes from pion cooling in the magnetic field, it can be computed from Eqs. (tsl) and (34) 



taking into account that about 1/4 of the pion energy is deposited per neutrino species. In 
order to reproduce Eiy^^^ ~ 10^ GeV, one has 5' ~ 3 ■ 10^ G. For the fiux normalization, we 
use Eq. ([I]) for this reference parameter set. We summarize these parameters in Table [I] 
(last column). 



However, as demonstrated in Ref. 14 , there are qualitative differences to this analytical 
picture because: 1. Multi-pion production processes in the photohadronic interactions affect 
the spectral shape. 2. Magnetic field effects, different for pions, muons, and kaons (c/.. 



Fig. 17), lead to a characteristic multi-peak spectrum. 3. Flavor mixing and pile- up effects 



pronounce the peak from muon decays. We will explicitely demonstrate these effects below. 

A quantitative example 

We show the steady pion (upper left) and muon (upper right) spectra, as well as the elec- 
tron neutrino (lower left), and muon neutrino (lower right) injection spectra for the chosen 
parameter values (c/.. Table III last column) in the SRF in Fig. 18 The steady pion spectra 



N'^[E'^) exhibit two breaks: one at about 10^ GeV, which comes from the break in the photon 
spectrum, and one at about 10^ GeV, determined by Eq. (34), above which the synchrotron 



losses dominate. Note that the pion injection and steady state spectra below the second 



break are connected with Eq. (35), which leads to a different spectral index by one power 
coming from t'^~cay; ^^^ -^^l- Y^\' "^^^ steady muon spectra N'^(E'^) are shown in the upper 



right panel of Fig. 18 Herepthere are in fact two cooling breaks visible at about 10 GeV 



(muons) and 10^ GeV (pions, from which the muons originate), see Fig. 17 Comparing 
the solid curves, including energy losses, with the dashed curves, without energy losses, the 
pile-up effect (the part of the spectrum where the flux including cooling is actually higher 
than without cooling) is more pronounced for the muons than for the pions. 



In the lower two panels of Fig. 18, we show the electron (left) and muon (right) neutrino 
injection spectra at the source, where the contributions from the individual parents are 
marked (c/.. Fig. [I]). Since electron neutrinos or antineutrinos (left panel) can only be 
produced in muon or neutron decays, the muon decays clearly lead to a pronounced peak. 
The neutron decays only contribute signiflcantly at very low energies. The muon neutrinos 
(right panel) can be produced from muon, pion, or kaon decays. One can easily see in this 
figure that the different critical energies in Eq. ^ (see also Fig. 17) lead to a characteristic 



triple-peak structure, where the kaons are least affected by the synchrotron cooling. 
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Figure 18: Pion (upper left), muon (upper right), electron neutrino (lower left), and niuon neutrino 
(lower right) spectra in the SRF, as given in the plot labels, for the reference parameters in Table IT] Solid 
curves include energy losses of the secondaries, dashed curves are shown without these energy losses. 
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Figure 19: Electron (left panel) and muon/tau (right panel) neutrino fluxes at the detector including 
Lorentz boost, redshift, and flavor mixing for the reference parameters in Table ^ (thick solid curves). The 
individual contributions from the spectra before flavor mixing are shown as well (thin solid curves). The 
neutrino energy flux density of the total fluxes (thick curves) is normalized to Eqs. M and (36), shown as 
"WB spectrum" . The shaded regions show the impact of the 3cr mixing angle uncertainties now, in about 
2015, and in about 2025 (see main text for details). The uncertainty in 2025 is basically smaller than the 
thickness of the curve. 



In Fig. [19} we include the effects of flavor mixing, Lorentz boost, and redsliift. Here, 
tlie electron (left panel) and muon or tau (right panel) neutrino fluxes at the detector 
(excluding Earth attenuation) for the reference parameters in Table fl] are shown (thick 
solid curves). The individual contributions from the spectra before flavor mixing are shown 
as well (thin solid curves), which correspond to the thick solid curves in the lower panels 
of Fig. flsl correctly weighted. The neutrino energy flux density of the total fluxes (thick 
curves) is normalized to Eqs. ([I]) and (36), shown as "WB spectrum". For the electron 
neutrino spectrum, a clear dominant peak can be seen, since the main contribution comes 
from the z/g spectrum before flavor mixing. For the muon/tau neutrino spectra, which are 
representative for muon tracks in neutrino telescopes, flavor mixing pronounces the peak 
from muon decays somewhat. However, the peak flux is dominated by a double peak at 
around the WB plateau, and the original shape is roughly reproduced. The additional kaon 
decay hump leads to a triple-peak spectrum. Note that even the neutrino spectrum from 
pion decays is not flat at the WB plateau, since high-energy photohadronic processes change 
the spectral shape, whereas the original shape can be exactly reproduced if these processes 
are switched off fT4|. 

Since in Eq. ^ the neutrino mixing angles appear, the uncertainties of the mixing angles 
will lead to uncertainties in the fluxes and flavor ratios. In order to discuss these, we use 
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Figure 20: Flavor ratio i? as a function of E^j at the detector. Here, "Pion beam" refers to a flavor ratio 
0e : 0p : (/)t = 1 : 2 : before flavor mixing, "Muon damped" to : 1 : 0, and "Neutron beam" to 1 : : 0. 
The dashed regions show the impact of the 3cr mixing angle uncertainties now ("current"), in about 2015, 
and in about 2025 (see main text for details). 



the following values and current Sa ranges (see Ref. [36]): sin 6*23 
"'^12 = 0.318 (3a: 0.27 ... 0.38), sin^ ^13 = (3a: sin^ ^13 < 



= 0.5 (3a: 0.36 ... 0.67), 
0.053). In addition, we 



sm' 

consider future improved bounds (3a) on ^13 from Daya Bay sin^ 6*13 < 0.012 and 6^23 from 
T2K 0.426 < sin^ ^23 < 0.574, which represent the expected sensitivities at around 2015 60 
if 6'i3 ~ 0, which we assume for the sake of simplicity. Beyond these, a neutrino factory may 
improve the bounds even further: sin^ ^13 < 1.5-10"^ (6lj and 0.46 < sin^^23 ^ 0.54 62 
(3a), which we label "2025". For Qyi-, the errors may improve somewhat, but new generations 
of experiments are not assumed here. Since Q\2 is of secondary importance in Eq. d6]), we 
do not consider such future improvements. The shaded regions in Fig. 19 show the ranges 
for the flux including these 3a uncertainties, where the different shadings correspond to the 
current errors, the ones in "2015", and the ones in "2025". Obviously, there is currently 
some uncertainty in the fluxes in the pion- and kaon-decay dominated ranges, where the 
v^ production dominates, especially for the v^ flux at the observer. However, already the 
improved knowledge from the upcoming experiments will drastically reduce this uncertainty. 
If information from a neutrino factory is available, the errors collapse to within the thick 
curves. 



Finally, we show in Fig. 20 the flavor ratio i? as a function of Ey at the detector. One can 



easily see the flavor transition from "pion beam", where the whole decay chain in Eq. (28) 



or Eq. (29) is present, to a "muon damped source", where the muons lose energy faster than 



they decay and the muon decay part is not present. In the presence of flavor and magnetic 
field effects, our results are consistent with Ref. 13 , in spite of the slightly different methods; 
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see also Ref. 33 . Again, we show the uncertainty coming from the mixing angles as dashed 



regions. Here, on a linear scale, the effect is somewhat larger than for the flux. However, 
again, in about 2015 the different regimes (pion beam and muon damped) can, in principle, 
be discriminated by the value of R. 

In Sec. [3| we use Fig. 19, right panel, and Fig. 20 as the representatives for the observables. 



A. 3 Limitations of the model, comparison to other approaches 



Our model can be regarded as minimal approach reproducing the results from Refs. 5 8,11 



which means that it also shares some of the weaknesses with these approaches. First of all, 
in all these approaches the proton spectrum and target photon field are put in "by hand" , 
which is motivated from observations, but may not be self-consistent for the individual 
burst. For example, the photohadronic cooling is not taken into account, compared to more 
detailed numerical simulations as in Ref. 12 . In the target photon spectrum, additional 



radiation processes may be at work, such as proton synchrotron radiation, which are not 



described by our approach (see, e.g.^ Refs. 47,63]). In addition, the effects of adiabatic 
cooling and additional interactions are neglected for the secondaries, which may change the 
spectral shape in extreme cases. 



Compared to Refs. 5 8,11 , our model treats the energy losses of the secondaries explicitely. 



which is important in the presence of strong magnetic fields to describe the spectrum and 
flavor ratios accurately, including pile-up effects which cannot be modeled in analytical 
treatments. Whereas additional processes affecting the primaries (protons, photons) may 
change the spectral shape, it is important to note that the flavor ratio relies on the descrip- 
tion of the secondaries, which mostly depends on (known) particle physics. In addition, 
we include additional pion production modes in the photohadronics, i.e., direct (t-channel) 
production, higher resonances, and multi-pion production, and the helicity-dependent muon 
decays, which have some impact on the flavor ratio. As it is explicitely demonstrated in 
Refs. 14,33 , it is fair to say that the A(1232) resonance (in the particle physics sense) 
hardly ever dominates the charged pion production in GRBs, see Fig. 5 in Ref. 33 . Al- 



though the methods are different, our ingredients are similar to Ref. 13 
able to reproduce and confirm, but perhaps far more efficient. 



which we were 



B Variation of luminosity and impact of gamma-ray detection 
threshold 



One may argue that the threshold for the 7-ray detection may lead to a selection bias in a 
stacking analysis. This is not taken into account in Sees. 3.1 and|4| which means that the 
low luminosity GRBs at high redshifts may be over-counted there. This threshold is (for 



Swift) shown as gray-shaded area in Fig. 1 of Kistler et al. 45 and roughly excludes the 



region Liso/{4:nd'j^{z)) < const in the I/iso--2-plane. Note that this detection threshold may 
potentially lead to a qualitative difference between the diffuse and quasi-diffuse (from the 
stacking using the gamma-ray counterpart) neutrino fluxes. 

We have tested the impact of this threshold for several assumptions of the luminosity distri- 
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Figure 21: In the left panel, we show the distribution of the number of observable GRBs as a function of 
redshift (histograms) with the luminosity being integrated out, where the gamma-ray detection threshold is 
taken into account. The two shown 10 000 burst samples differ in their distribution in luminosity, with the 
dark one being wider (cr^j^ = 1.0) than the light one being narrower (ax2 — 0.5), see main text. The dashed 
curves depict the analytical result of the distribution in redshift, the thick solid curves the corresponding 
contribution functions (red/dark: wide, orange/light: narrow). Notice that the contribution functions are 
in a.u. and have been normalized to 1 when integrated over the depicted redshift range. In the right 
panel, the cumulative distribution function similar to Fig. [15] is shown for n = 100 and n = 10 000 for 
only redshift varied (dashed) and redshift and luminosity varied simultaneously (solid), where the wide 
luminosity distribution is used. 



bution function by simultaneously generating redshift and luminosity distributions. First of 
all, it should be noted that we do not observe a significant impact of the threshold on the con- 
tribution function or neutrino spectra as long as the stacked bursts satisfy Liso ^ 10^^ erg s~\ 
i.e., a cut on the luminosity is applied. Since this cut is exactly the one which has been 
also used in Ref. Es] to derive the correction to the star formation rate evolution of the 
Swift GRBs, our redshift distribution is based on, this result is not surprising. Note that, of 
course, the cut on the luminosity depends on the satellite used for the stacking sample. As 
a qualitative difference, if an additional luminosity distribution is applied, the cumulative 
distribution function of the bursts will be broader (the statistical errors add, to a first ap- 
proximation, in a Gaussian way), which means that the systematical errors given in Sec. 
on the stacked flux can be regarded as a lower limit. The details, however, depend on the 
assumed luminosity distribution, see below for an example, which is why we do not include 
the luminosity in the main text. 



If bursts with Liso ^ 



10 ergs" are included in the stacking analysis, i.e., bursts which 
may be visible for small z but not for high z, the luminosity distribution may have an 
impact. Assume, for instance, that the GRB luminosities follow a Gaussian distribution in 
the exponent Liso = lO^ergs"^ and a central value x ~ 51, see Fig. 1 in Ref. 45 . Consider 
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a wider standard deviation ax^ = 1-0, and a narrower one ax2 = 0.5 for comparisonF^ Both 
of these functions are expected to have a significant effect on the result of the observed 
number of bursts compared to Fig. |3] if the Swift detection threshold is applied. We show 



the results in Fig. [21} In the left plot, the redshift distributions with luminosity integrated 
out are shown, together with the corresponding contribution functions. As one can see from 
the plot, the high redshift contribution is suppressed for the narrower distribution, because 
not as many bursts exceeding the threshold as for the wide distribution are generated 
in that region. However, the luminosity-integrated contribution functions both peak at 
2; ~ 1 and are very similar to the one in Fig. |3j As far as the cumulative distribution 



is concerned (right panel of Fig. 21 for the wide distribution function of luminosity), the 
spreads including the luminosity distribution (solid curves) are significantly larger than for 
the redshift distribution (dashed curves) only, which is expected from statistics. In fact, in 
this example, the 90%CL error on the stacking of 100 bursts is roughly 100%. That number, 
however, depends on the breadth of the luminosity distribution of the stacked bursts, which 
is chosen aggressively wide in this example. In either case, the redshift distribution will 
limit this systematical error, as we argue in the main text. 



-•^^There have been some studies to extract the luminosity distribution of the Swift sample, see Refs. 64|65 



For instance, the broken power law distribution in Ref. i64l leads to results between our two extreme cases. 

45 



